cis-meQTLs analysis: Genotype with missingness. LD clumped results. meQTL vs. eQTM

Dexamethasone-Stimulated Human Array Project

Max Planck Institute of Psychiatry
anastasiia_hry@psych.mpg.de

Introduction

meQTLs

This markdown provides comparison between two meQTL analysis after LD clumping. meQTL analysis were done using the R package MatrixEQTL. LD clumping was done using PLINK

For all analysis, the same DNAm data were used which underwent following QC steps:

DNAm QC

The imputed QC SNP data (no LD pruning, the MHC regions are included) were used. The data includes 3’957’338 SNPs.

All meQTLs passed the significance threshold of FDR = 0.05.

LD clumping was done per CpG using the following parameters:

–clump-p1 = 0.05: Significance threshold for index SNPs

–clump-p2 = 1: Secondary significance threshold for clumped SNPs

–clump-r2 = 0.2: LD threshold for clumping

–clump-kb = 200: Physical distance threshold for clumping

eQTMs

Results

Overlap meQTLs & eQTMs

Baseline

Overlap between baseline meQTL CpGs (pre-dexamethasone) and baseline eQTM CpGs (pre-dexamethasone)

euler.tbl <- euler(euler.cpgs)

# For visualisation purposes:
euler.tbl$ellipses["meqtl", c("a", "b")] <- euler.tbl$ellipses["meqtl", c("a", "b")]

plot(euler.tbl, quantities = list(type = c("counts", "percent")), fills = list(fill = c("grey", "transparent"),
    alpha = 0.5), labels = c("meQTL", "eQTM"))

intersect.cpgs <- intersect(euler.cpgs$meqtl, euler.cpgs$eqtm)

eqtm.neg.cor <- eqtm.veh.nom.df[CpG_ID %in% intersect.cpgs][beta_eqtm < 0] %>%
    nrow()
eqtm.pos.cor <- eqtm.veh.nom.df[CpG_ID %in% intersect.cpgs][beta_eqtm > 0] %>%
    nrow()

perc.eqtm.neg.cor <- percent(eqtm.neg.cor/nrow(eqtm.veh.nom.df[CpG_ID %in% intersect.cpgs]), accuracy = 0.01)
perc.eqtm.pos.cor <- percent(eqtm.pos.cor/nrow(eqtm.veh.nom.df[CpG_ID %in% intersect.cpgs]), accuracy = 0.01)

• 42% (57675) meQTL CpGs overlap with baseline (pre-dexamethasone) eQTM CpGs: 52.06% are negatively and 47.94% positively correlated.

• 57675 meQTL CpGs are mapped to 11381 eQTM ENSG (152763 eQTMs).

GR-induced

Overlap between GR-induced meQTL CpGs (post-dexamethasone - pre-dexamethasone) and GR-induced eQTM CpGs (post-dexamethasone)

euler.tbl <- euler(euler.cpgs)

# For visualisation purposes:
euler.tbl$ellipses["meqtl", c("a", "b")] <- euler.tbl$ellipses["meqtl", c("a", "b")] * 5

plot(euler.tbl, quantities = list(type = c("counts", "percent")), fills = list(fill = c("#0072B2", "transparent"),
    alpha = 0.5), labels = c("meQTL", "eQTM"))

intersect.cpgs <- intersect(euler.cpgs$meqtl, euler.cpgs$eqtm)

eqtm.neg.cor <- eqtm.veh.nom.df[CpG_ID %in% intersect.cpgs][beta_eqtm < 0] %>%
    nrow()
eqtm.pos.cor <- eqtm.veh.nom.df[CpG_ID %in% intersect.cpgs][beta_eqtm > 0] %>%
    nrow()

perc.eqtm.neg.cor <- percent(eqtm.neg.cor/nrow(eqtm.veh.nom.df[CpG_ID %in% intersect.cpgs]), accuracy = 0.01)
perc.eqtm.pos.cor <- percent(eqtm.pos.cor/nrow(eqtm.veh.nom.df[CpG_ID %in% intersect.cpgs]), accuracy = 0.01)

• 1% (903) meQTL CpGs overlap with GR-induced (post-dexamethasone) eQTM CpGs: 55.13% are negatively and 44.87% positively correlated.

• 903 meQTL CpGs are mapped to 1352 eQTM ENSG (1785 eQTMs).

Heatmap

corrplot(meth.gex.cor.mtrx, method = "color", cl.pos = "r", type = "full", tl.pos = "n", title = paste0("Top ",
    n.top.eqtm, " eQTMs, ", nrow(meth.sub.df), "CpG,  ", nrow(gex.sub.df), "GEX"))

Correlation within cis-window of 1Mbp, CpG and GEX delta with parallel FC

Examples

Only post-dex

intersect.cpgs <- intersect(meqtl.delta.df$CpG_ID, meqtl.veh.df$CpG_ID)
selected.meqtl <- meqtl.delta.df[!(CpG_ID %in% intersect.cpgs)][fdr == min(fdr)]

selected.eqtm <- eqtm.dex.nom.df[CpG_ID %in% selected.meqtl$CpG_ID, .(CpG_ID, ENSG_ID, beta_eqtm, `p-value_eqtm`,
    fdr_eqtm)] %>%
    unique()

selected.eqtl <- data.frame(ENSG_ID = selected.eqtm$ENSG_ID, SNP = selected.meqtl$SNP) %>%
    unique()

meQTL

data.table(selected.meqtl %>%
    dplyr::select(CpG_ID, SNP, FC = beta, p_val = `p-value`, FDR = fdr) %>%
    mutate_if(is.numeric, funs(as.character(signif(., 3)))))

ProcessGetBoxPlot(meth.beta.veh.mtrx, meth.beta.dex.mtrx, snp.mtrx, selected.meqtl[1], fdr.thr = 0.05,
    plot.title = plot.title)

eQTM

data.table(selected.eqtm %>%
    dplyr::select(CpG_ID, ENSG_ID, FC = beta_eqtm, p_val = `p-value_eqtm`, FDR = fdr_eqtm) %>%
    mutate_if(is.numeric, funs(as.character(signif(., 3)))))

ScatterPlotGEXvsDNAm(meth.beta.dex.mtrx, meth.beta.veh.mtrx, gex.dex.mtrx, gex.veh.mtrx, cpg.id = selected.eqtm$CpG_ID[1],
    ensg.id = selected.eqtm$ENSG_ID[1])

eQTL

ProcessGetBoxPlot(gex.veh.mtrx, gex.dex.mtrx, snp.mtrx, selected.eqtl[1, ], fdr.thr = 0.05, plot.title = plot.title)

Both pre and post-dex

intersect.cpgs <- intersect(meqtl.delta.df$CpG_ID, meqtl.veh.df$CpG_ID)
selected.meqtl <- meqtl.delta.df[CpG_ID %in% intersect.cpgs][fdr == min(fdr)]

selected.eqtm <- eqtm.dex.nom.df[CpG_ID %in% selected.meqtl$CpG_ID, .(CpG_ID, ENSG_ID, beta_eqtm, `p-value_eqtm`,
    fdr_eqtm)] %>%
    unique()

selected.eqtl <- data.frame(ENSG_ID = selected.eqtm$ENSG_ID, SNP = selected.meqtl$SNP) %>%
    unique()

meQTL

data.table(selected.meqtl %>%
    dplyr::select(CpG_ID, SNP, FC = beta, p_val = `p-value`, FDR = fdr) %>%
    mutate_if(is.numeric, funs(as.character(signif(., 3)))))

ProcessGetBoxPlot(meth.beta.veh.mtrx, meth.beta.dex.mtrx, snp.mtrx, selected.meqtl[1], fdr.thr = 0.05,
    plot.title = plot.title)

eQTM

data.table(selected.eqtm %>%
    dplyr::select(CpG_ID, ENSG_ID, FC = beta_eqtm, p_val = `p-value_eqtm`, FDR = fdr_eqtm) %>%
    mutate_if(is.numeric, funs(as.character(signif(., 3)))))

ScatterPlotGEXvsDNAm(meth.beta.dex.mtrx, meth.beta.veh.mtrx, gex.dex.mtrx, gex.veh.mtrx, cpg.id = selected.eqtm$CpG_ID[1],
    ensg.id = selected.eqtm$ENSG_ID[1])

eQTL

ProcessGetBoxPlot(gex.veh.mtrx, gex.dex.mtrx, snp.mtrx, selected.eqtl[1, ], fdr.thr = 0.05, plot.title = plot.title)